{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f5910b3b",
   "metadata": {},
   "source": [
    "这是一个关于零极点分析的代码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "31ef59ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUXklEQVR4nO3df7DddX3n8edLIu2UqMCikV82uEu7ws7YSoK4a0vu+GOB1knrIMXtKu6skwHL7o6zDsW12Fb3h2VmO1vLL2nHFtapMZ3VmtIoq+wN1JnSDanGEpCSpVDS8KMgPwzW2uB7/zjf1JvLvTfn5pyce7/5PB8zd+731/l+XvnCvM73fM/3npOqQpJ05HvRUgeQJE2GhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX8tOkvck+coC67+Q5JJJZhqXJDckueow7Le3x0STY+FrIpKsTrIlyVNJHk1yTZIVh7Kvqjq/qm7q9rvgk8NSmitbVV1aVR8dcb+/kuRTs/b7D8dEmo+Fr0m5DngcOBH4MeBc4H1LGUhqjYWvSTkN2FRV36mqR4EvAmcusH2S/GaSZ5J8I8mbZqzYmuS9SV4D3AC8IcneJE/Ps6PTktye5FtJvtS9uvhUt25dkt2ztn8wyZu76bOT/EmSp5M80j326BnbVpJLk9zfvXq5NgNzZkvyu0n+czf9h926/T/fS/Kebt1vJHk4ybNJtif5iW75ecB/An6ue8yOmcekm35Rkl9K8lCSx5PcnORl3brVXeZLkvxVkieSfGiI/346Alj4mpTfAC5O8kNJTgbOZ1D683k98ABwAvDLwGeTHD9zg6q6F7gU+JOqWllVx86zr98Dtnf7+iiwmGvdzwPv7x77BuBNvPCVyU8Da4HXAhcB/3KYbFX1tm7dSuBC4FHgtm71NgavhI7v8v9+kh+sqi8C/xX4TPfY186R+T3dzxTwamAlcM2sbd4I/Gj37/lw9wSlI5yFr0m5ncEZ/bPAbuAu4A8W2P5x4H9U1d9X1WeA+4CfWuygSV7FoIyvqqq/q6o7gD8c9vFVtb2q7qyqfVX1IPAJBpejZvpYVT1dVX8FTDMo6sVk/BHgZuDnqurhbtxPVdWT3bj/HfgBBgU9jJ8Hfr2qHqiqvcAHGTzZznzP5Fer6m+ragewg8GTlY5wFr4OuyQvAm4FPgscw+Bs+Tjg1xZ42F/XgZ/s9xBw0iEMfxLwVFU9N2tfQ0nyI0lu6d5ofpbB2fUJszZ7dMb0txmcUQ+7/5cBn2fwhPTHM5b/xyT3dpe0ngZeNse48zmJA/+NDwErgFXjyKz+svA1CccDpwLXdGfZTwK/A1ywwGNOTpIZ868C9syx3cE+7vUR4Lgkx8za137PAT+0fybJUcDLZ6y/HvgGcHpVvZTB9fOZuRayYLbuifD3gOmq+sSM5T8B/CKDy0PHdZeDnpkx7sH+zXuAH54x/ypgH/DYkLl1hLLwddhV1RPAXwKXJVmR5FgG19F3LPCwVwD/PsmLk7wDeA2wZY7tHgNOmflG6qyxH2Jw+ehXkxyd5I3A22Zs8hfADyb5qSQvBn6JweWT/V7C4DLU3iT/FLjs4P/i4bIB/4XBK57/MGv5SxgU9N8AK5J8GHjprP2u7p4w5vJp4P3dm9Ur+f41/32LyK4jkIWvSXk7cB6DEtvFoNDev8D2fwqcDjzBoBgv7F4ZzPZ/gJ3Ao0memGdf/4rBm8DfZPAG8M37V1TVMwzehP1t4K8ZnPHPvGvnA93jvwX8FvCZhf6Ri8z2TuAc4KkZd+r8PIPLX19g8GT0EPAd4OEZj/v97veTSf5sjv1+EvifwB0Mnmi/A/y7ReTWESp+AYpak+RXgH9SVf96qbNIk+QZviQ1wsKXpEZ4SUeSGuEZviQ14pA+rXBSTjjhhFq9evVExnruuec45phjDr7hMtTX7H3NDf3N3tfc0N/sk869ffv2J6rq5XOtW9aFv3r1au66666JjLV161bWrVs3kbHGra/Z+5ob+pu9r7mhv9knnTvJvH9J7iUdSWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRoyl8JN8MsnjSe6eZ32SfDzJriRfT/K6cYwrSRreuM7wfxc4b4H15wOndz8bgOvHNK4kaUhjKfyqugP45gKbrAduroE7gWOTnDiOsSVJwxnbd9omWQ3cUlX/bI51twAfq6qvdPO3Ab9YVS/4dpMkGxi8CmDVqlVnbdy4cSz5Dmbv3r2sXLlyImONW1+z9zU39Dd7X3NDf7NPOvfU1NT2qloz17pJfeNV5lg25zNNVd0I3AiwZs2amtQ3xfT123Sgv9n7mhv6m72vuaG/2ZdT7kndpbMbOHXG/CnAngmNLUlicoW/GXh3d7fOOcAzVfXIhMaWJDGmSzpJPg2sA05Ishv4ZeDFAFV1A7AFuADYBXwb+DfjGFeSNLyxFH5VvfMg6wv4hXGMJUk6NP6lrSQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiPGUvhJzktyX5JdSa6cY/26JM8k+Vr38+FxjCtJGt6KUXeQ5CjgWuAtwG5gW5LNVXXPrE3/uKp+etTxJEmHZhxn+GcDu6rqgar6LrARWD+G/UqSxihVNdoOkguB86rqvd38u4DXV9XlM7ZZB/wvBq8A9gAfqKqd8+xvA7ABYNWqVWdt3LhxpHzD2rt3LytXrpzIWOPW1+x9zQ39zd7X3NDf7JPOPTU1tb2q1sy5sqpG+gHeAfz2jPl3Ab85a5uXAiu76QuA+4fZ91lnnVWTMj09PbGxxq2v2fuau6q/2fuau6q/2SedG7ir5unUcVzS2Q2cOmP+FAZn8TOfVJ6tqr3d9BbgxUlOGMPYkqQhjaPwtwGnJzktydHAxcDmmRskeWWSdNNnd+M+OYaxJUlDGvkunaral+Ry4FbgKOCTVbUzyaXd+huAC4HLkuwD/ha4uHvpIUmakJELH/7hMs2WWctumDF9DXDNOMaSJB0a/9JWkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsJvzNVXw/T0gcumpwfLJR3ZLPzGrF0LF130/dKfnh7Mr127tLkkHX4rljqAJmtqCjZtGpT8ZZfB9dcP5qemljqZpMPNM/wGTU0Nyv6jHx38tuylNlj4DZqeHpzZX3XV4Pfsa/qSjkwWfmP2X7PftAk+8pHvX96x9KUjn4Xfkquv5rGN0wdcs59imjvffjXbti1tNEmHn4XfkrVrufizFzHFgbfo/OOL13LFFUsbTdLh5106LfEWHalpnuG3xlt0pGZZ+K3xFh2pWRZ+S7xFR2qahd+SbdsOvGa//5q+t+hITfBN25bMdSvO1JTX8aVGjOUMP8l5Se5LsivJlXOsT5KPd+u/nuR14xhXmrStZ76PfVnBuVNT7MsKtp75vqWOJA1t5MJPchRwLXA+cAbwziRnzNrsfOD07mcDcP2o40qTtvXM93HuPdezgucJsILnOfee6y199cY4zvDPBnZV1QNV9V1gI7B+1jbrgZtr4E7g2CQnjmFsaWLeeM+NZNaydMulPhjHNfyTgYdnzO8GXj/ENicDj8zeWZINDF4FsGrVKrZu3TqGiAe3d+/eiY01bn3N3rfc5/L8nMuP4vne/Dv6dsxn6mv25ZR7HIU/+6QHoA5hm8HCqhuBGwHWrFlT69atGyncsLZu3cqkxhq3vmbvW+59HMWKOUr/eY7qzb+jb8d8pr5mX065x3FJZzdw6oz5U4A9h7CNtKx95YwNLzhLqW651AfjKPxtwOlJTktyNHAxsHnWNpuBd3d365wDPFNVL7icIy1n63Zex+1nXMY+jqIYnPHffsZlrNt53VJHk4YycuFX1T7gcuBW4F5gU1XtTHJpkku7zbYADwC7gN8CvK1BvbRu53WsqH3cPj3Nitpn2atXxvKHV1W1hUGpz1x2w4zpAn5hHGNJkg6NH63QkquvfuHn5kxPD5ZLOuJZ+C1Zu/bAD0vb/2Fqa9cubS5JE+Fn6bTEL0CRmuYZfmv8AhSpWRZ+a/wCFKlZFn5Lpqf59tsu4qsfPPALUL7669O+bys1wMJvybZt3PeRTbz1v00NTuynpvjqBzex+aptvm8rNcA3bVtyxRX8OLDpx2e+bzvFplumvJQvNcAz/Ab5vq3UJgu/Qb5vK7XJwm/M/r+12nTg+7aWvtQAC78x27Yd+LdW+/8Wa9u2pc0l6fDzTdvGXHHFC5dNTXkdX2qBZ/iS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNWLFKA9OcjzwGWA18CBwUVU9Ncd2DwLfAp4H9lXVmlHGlSQt3qhn+FcCt1XV6cBt3fx8pqrqxyx7SVoaoxb+euCmbvom4GdG3J8k6TBJVR36g5Onq+rYGfNPVdVxc2z3l8BTQAGfqKobF9jnBmADwKpVq87auHHjIedbjL1797Jy5cqJjDVufc3e19zQ3+x9zQ39zT7p3FNTU9vnvZJSVQv+AF8G7p7jZz3w9Kxtn5pnHyd1v18B7AB+8mDjVhVnnXVWTcr09PTExhq3vmbva+6q/mbva+6q/mafdG7grpqnUw/6pm1VvXm+dUkeS3JiVT2S5ETg8Xn2saf7/XiSzwFnA3ccbGxJ0viMeg1/M3BJN30J8PnZGyQ5JslL9k8Db2XwCkGSNEGjFv7HgLckuR94SzdPkpOSbOm2WQV8JckO4P8Cf1RVXxxxXEnSIo10H35VPQm8aY7le4ALuukHgNeOMo4kaXT+pa0kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjRir8JO9IsjPJ95KsWWC785Lcl2RXkitHGVOSdGhGPcO/G3g7cMd8GyQ5CrgWOB84A3hnkjNGHFeStEgrRnlwVd0LkGShzc4GdlXVA922G4H1wD2jjC1JWpxU1eg7SbYCH6iqu+ZYdyFwXlW9t5t/F/D6qrp8nn1tADYArFq16qyNGzeOnG8Ye/fuZeXKlRMZa9z6mr2vuaG/2fuaG/qbfdK5p6amtlfVnJfYD3qGn+TLwCvnWPWhqvr8EOPPdfo/77NMVd0I3AiwZs2aWrdu3RBDjG7r1q1Maqxx62v2vuaG/mbva27ob/bllPughV9Vbx5xjN3AqTPmTwH2jLhPSdIiTeK2zG3A6UlOS3I0cDGweQLjSpJmGPW2zJ9Nsht4A/BHSW7tlp+UZAtAVe0DLgduBe4FNlXVztFiS5IWa9S7dD4HfG6O5XuAC2bMbwG2jDKWJGk0/qWtJDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUCAtfkhoxlm+8OlyS/A3w0ISGOwF4YkJjjVtfs/c1N/Q3e19zQ3+zTzr3D1fVy+dasawLf5KS3DXf14Itd33N3tfc0N/sfc0N/c2+nHJ7SUeSGmHhS1IjLPzvu3GpA4ygr9n7mhv6m72vuaG/2ZdNbq/hS1IjPMOXpEZY+JLUiGYLP8k7kuxM8r0k894yleTBJH+e5GtJ7ppkxvksIvt5Se5LsivJlZPMOE+e45N8Kcn93e/j5tlu2Rzzgx3DDHy8W//1JK9bipyzDZF7XZJnumP8tSQfXoqcsyX5ZJLHk9w9z/plebxhqOxLf8yrqskf4DXAjwJbgTULbPcgcMJS511sduAo4P8BrwaOBnYAZyxx7quBK7vpK4FfW87HfJhjCFwAfAEIcA7wpz3JvQ64ZamzzpH9J4HXAXfPs37ZHe9FZF/yY97sGX5V3VtV9y11jkMxZPazgV1V9UBVfRfYCKw//OkWtB64qZu+CfiZpYsylGGO4Xrg5hq4Ezg2yYmTDjrLcvxvP5SqugP45gKbLMfjDQyVfck1W/iLUMD/TrI9yYalDrMIJwMPz5jf3S1bSquq6hGA7vcr5tluuRzzYY7hcjzOw2Z6Q5IdSb6Q5MzJRBvZcjzei7Gkx3zFpAecpCRfBl45x6oPVdXnh9zNv6iqPUleAXwpyTe6Z/LDagzZM8eyw34P7kK5F7GbJTnmcxjmGC7JcT6IYTL9GYPPXNmb5ALgD4DTD3ewMViOx3tYS37Mj+jCr6o3j2Efe7rfjyf5HIOXy4e9fMaQfTdw6oz5U4A9I+7zoBbKneSxJCdW1SPdy/DH59nHkhzzOQxzDJfkOB/EQTNV1bMzprckuS7JCVW13D+cbDke76Esh2PuJZ0FJDkmyUv2TwNvBeZ8B34Z2gacnuS0JEcDFwOblzjTZuCSbvoS4AWvVJbZMR/mGG4G3t3dPXIO8Mz+y1ZL6KC5k7wySbrpsxl0wZMTT7p4y/F4D2VZHPOlfmd7qX6An2VwtvB3wGPArd3yk4At3fSrGdzhsAPYyeBySi+yd/MXAH/B4I6NJc8O/CPgNuD+7vfxy/2Yz3UMgUuBS7vpANd26/+cBe74Wma5L++O7w7gTuCfL3XmLtengUeAv+/+H/+3fTjeQ2Zf8mPuRytIUiO8pCNJjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiP+P//fIBUHSWt9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from scipy import signal\n",
    "\n",
    "from matplotlib.patches import Circle\n",
    "\n",
    " \n",
    "N = 8 # quantization bits\n",
    "\n",
    "Fs = 1000 # sampling frequency\n",
    "\n",
    "a = np.array([1, 1.7, 0.745]) # denominator\n",
    "\n",
    "b = np.array([0.05, 0, 0]) # numerator\n",
    "\n",
    " \n",
    "# z zero, p pole, k gain\n",
    "\n",
    "z1, p1, k1 = signal.tf2zpk(b,a) # zero, pole and gain\n",
    "\n",
    " \n",
    "c = np.vstack((a,b))\n",
    "\n",
    "Max = (abs(c)).max() #find the largest value\n",
    "\n",
    "a = a / Max #normalization\n",
    "\n",
    "b = b / Max\n",
    "\n",
    " \n",
    "Ra = (a * (2**((N-1)-1))).astype(int) # quantizan and truncate\n",
    "\n",
    "Rb = (b * (2**((N-1)-1))).astype(int)\n",
    "\n",
    "z2, p2, k2 = signal.tf2zpk(Rb,Ra)\n",
    "\n",
    " \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "circle = Circle(xy = (0.0, 0.0), radius = 1, alpha = 0.9, facecolor = 'white')\n",
    "\n",
    "ax.add_patch(circle)\n",
    "\n",
    "for i in p1:\n",
    "\n",
    " ax.plot(np.real(i),np.imag(i), 'bx') #pole before quantization\n",
    "\n",
    "for i in z1:\n",
    "\n",
    " ax.plot(np.real(i),np.imag(i),'bo') #zero before quantization\n",
    "\n",
    "for i in p2:\n",
    "\n",
    " ax.plot(np.real(i),np.imag(i), 'rx') #pole after quantization\n",
    "\n",
    "for i in z2:\n",
    "\n",
    " ax.plot(np.real(i),np.imag(i),'ro') #zero after quantization\n",
    "\n",
    " \n",
    " \n",
    "ax.set_xlim(-1.8,1.8)\n",
    "\n",
    "ax.set_ylim(-1.2,1.2)\n",
    "\n",
    "ax.grid()\n",
    "\n",
    "ax.set_title(\"%d bit quantization\" %N)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
